#ifndef __events
#define __events
#include"particles.h"
#include "physicalconstants.h"
#include"type.h"
#include "userdef.h"
#include <random>
#include <vector>
#include <algorithm>
///说明：events调用具体的mc过程，然后对新粒子的od初始化和随机化
/// 在粒子中封装events，具体类型的sp调用不同的过程函数
/// 不能使用map，因为采用map时虽然产找方便，但变量的取值就出问题了。
/// 这个时外部接口,通过trigger,可以从这个对象中获得新的粒子，以及修改过的粒子

// 应该吧这两个函数放到common中。
/// return the index of param in matrix
//template<typename __Ty> 
//int bs_search(const __Ty& param, const vector<__Ty>& matrix)
int bs_search(double& param, vector<double>& matrix)
{
    vector<double>::iterator locator;
    locator = lower_bound(matrix.begin(), matrix.end(), param);
    int ret = (locator - matrix.begin()) - 1;
    if(ret == matrix.size())
    {
        ret = - 1;
    }
    if(ret == 0)
    {
        ret = -1;
    }
    return ret;
}

double init_velocity(const double& lorentz_gamma)
{
    return c_light_speed * sqrt(1.0 - 1.0 / lorentz_gamma / lorentz_gamma);
}

double random_generator()
{
  std::random_device rd;
  std::mt19937 gen(rd());
  std::uniform_real_distribution<> distribution(0.0, 1.0);
  return distribution(gen);
}


void init_ptcls(__Ptcls& myptcl)
{
  myptcl.od_factor = 0.0;
  myptcl.rod_factor = - log(1.0 - random_generator());
}

/// 得到该矢量方向的单位向量
template<typename __Ty>
__Vect3<__Ty> normalize(__Vect3<__Ty>& normalizer)
{
    return normalizer / Mod(normalizer);
}

///更新正负电子和光子的QED参数
void update_qed_factor(__Ptcls& particle)
{
    switch(particle.species)
    {
        case 0:
        case 1:
            {
                particle.qed_factor = particle.lorentz_gamma * Mod(particle.efield - \
                        (particle.velocity * particle.efield) * particle.velocity / \
                        (particle.velocity * particle.velocity) + Cross(particle.velocity,\
                            particle.bfield)) / E_crit;
            }
            break;
        case 3:
                particle.qed_factor = particle.energy * Mod(particle.efield - \
                        (particle.velocity * particle.efield) * particle.velocity / \
                        (particle.velocity * particle.velocity) + Cross(particle.velocity, \
                                particle.bfield)) / E_crit / 2.0 / mc2;
            break;
        default:
            break;
    }
}



// 可能有2个表0,eta, 1, chi, 2 f, 3 h(eta), 4 T(chi), 5 chi for T, 6 eta for h
// 0 Pchi(eta, chi)  1, Pf(f, chi), 2: chi(eta, P)
// 返回值代表新的粒子个数
int synchrotron_rad(const double& dt, __Ptcls& parent, __Ptcls& child,\
        vector<vector<double> >& table, vector<vector<vector<double> > >& table_P)
{
    //cout<<"in syn"<<endl;
    //cout<<"parent beta is "<<Mod(parent.velocity) / c_light_speed<<" field "\
        <<parent.efield<<endl;
    double local_const = 9.81509e+18;
    if(parent.od_factor <= 1e-20)
    {
        //cout<<"initing rod"<<endl;
        init_ptcls(parent);
    }
    update_qed_factor(parent);
    /// 对于电子和正电子qed factor是ｅｔａ, 光子的是chi
    double temp_qed = log10(parent.qed_factor);
    /// 先得到eta的位置，利用eta的位置得到P，eat的位置处的chi
    int eta_index = bs_search(temp_qed, table[0]);
    int eta_index_h = bs_search(temp_qed, table[6]);
    
    if(eta_index < 0 || eta_index > table[0].size())
    {
        return 0;
    }
    
    //cout<<"eta index "<<eta_index_h<<" heta "<<table[3][eta_index_h]<<endl;
    parent.od_factor += dt * pow(10.0, table[3][eta_index_h]) * parent.qed_factor / parent.lorentz_gamma * local_const;
    //cout<<"local const is "<<sqrt(3.0) * alpha_fine * c_light_speed / lambda_compton;
    
    //cout<<"od factor now "<<parent.od_factor<<" rod is "<<parent.rod_factor<<endl;
    if(parent.od_factor < parent.rod_factor)
    {
        return 0;
    }
    int max_indP = table_P[0][eta_index].size() - 1;
    double random_P = random_generator() * table_P[0][eta_index][max_indP];
    // Ｐ是chi 和eta的函数
    int chi_index = bs_search(random_P, table_P[0][eta_index]);
    if( chi_index < 0 || chi_index > table_P[2][eta_index].size() - 1)
    {
        return 0;
    }
    double chi_f = pow(10.0, table_P[2][eta_index][chi_index]);
    /// 通过从父粒子获取信息，因为很多东西共享
    //很明显ＩＤ这东西没法通用了
    child = parent;
    child.species = 3;
    child.charge = 0.0;
    child.mass = 0.0;
    //cout<<" parent "<<parent<<endl;
    double pold_e = parent.energy;
    child.energy = 2.0 * mc2 * chi_f * parent.lorentz_gamma / parent.qed_factor;
    if(child.energy > parent.energy)
    {
        return 0;
    }
    //cout<<"parent "<<parent.species<<"\t"<<parent.energy/MeV<<\
        " MeV_____photon generated with ek "\
        <<child.energy/MeV<<" MeV______"<<"\t"<<parent.lorentz_gamma<<endl;
    child.velocity = c_light_speed * normalize(parent.velocity);
    update_qed_factor(child);
    init_ptcls(child);
    init_ptcls(parent);
    /// 这个地方计算正负电子的速度应该根据动量计算，一会儿再修改：
    /// p = (1 - hv / cp ) P
    //parent.energy = parent.energy - child.energy;
    //parent.lorentz_gamma = parent.energy / mc2;
    //cout<<" new velocity ratio is "<<(1.0 - child.energy / \
            (c_light_speed * Mod(parent.velocity)))<<endl;
    double old_gamma = parent.lorentz_gamma;
    parent.energy = parent.energy - child.energy;
    parent.lorentz_gamma = parent.energy / mc2 + 1.0;
    parent.velocity = parent.velocity * old_gamma / parent.lorentz_gamma * (1.0 - child.energy / old_gamma / \
            parent.mass / c_light_speed / Mod(parent.velocity));
    /// 需要设定光子能量的一个下限
    if(parent.energy < 0)
    {
        cout<<" parent negative energy "<<parent.energy<<endl;
        cout<<" parent e "<<parent.efield<<" parent v "<<parent.velocity<<endl;
        cout<<"parent qed "<<parent.qed_factor<<" chif "<<chi_f<<" parent old energy "<<pold_e<<" old gamma "<<old_gamma<<" c ek"<<child.energy<<endl;
        getchar();
    }
    if(child.energy < ekmin)
    {
        return 0;
    }
    if(std::isnan(parent.energy))
    {
        cout<<parent<<endl;
        cout<<" parent nan energy "<<endl;
    }
    return 1;
}





/// tables1d 0: ek_ele,1: x = omega / ek
/// tables2d: 0 sigma_br(ek_ele, x), 1, sigma_bh(ek_photon, y)
int bremsstrahlung_rad(const double& dt, const double& density, \
        __Ptcls& parent, __Ptcls& child, vector<vector<double> >& tables1d, \
        vector<vector<vector<double> > >& tables2d)
{
    
    int energy_index = bs_search(parent.energy, tables1d[0]);
    //cout<<"ek index "<<energy_index<<"\t"<<parent.energy<<"\t"<<parent.id<<"\t"<<parent.species<<endl;
    //getchar();
    if(energy_index < 0)
    {
        return 0;
    }
    /// offgrid 能量，采用区间端点概率决定用那个e
    //energy_index = (random_generator() > (log(tables1d[0][energy_index + 1]) - \
                log(parent.energy)) / (log(tables1d[0][energy_index + 1]) - \
                log(tables1d[0][energy_index])))?energy_index:(energy_index + 1);
    //cout<<" ph f_index "<<energy_index<<endl;
    int max_ind = tables1d[0].size();
    double sigma_total = tables2d[0][energy_index][max_ind - 1];
    /// cross section units is cm2
    double freq = Mod(parent.velocity) * density * sigma_total * dt;
    double probability = 1.0 - exp(-freq);
    //cout<<" parent ek "<<parent.energy<<" in table "<<tables1d[0][energy_index]<<endl;
    //cout<<"prob is "<<probability<<" with freq "<<freq<<" density "<<density<<"vel"<<Mod(parent.velocity)<<endl;
    double myrand = random_generator();
    if(myrand > probability)
    {
        return 0;
    }
    //cout<<"prob is "<<probability<<" myrand is "<<myrand<<endl;
    //cout<<"may be will emit photons"<<endl;
    double sigma = sigma_total * random_generator();
    int ph_index = 0; 
    while(ph_index <= 0)
    {
        sigma = sigma_total * random_generator();
        ph_index = bs_search(sigma, tables2d[0][energy_index]);
    }
    //cout<<"ph_index is "<<ph_index<<" ek f_index "<<energy_index<<endl;
    child = parent;
    /// 对ph_index不在格点的sigma进行插值。
    double ek_ratio = tables1d[1][ph_index] + \
                      (tables1d[1][ph_index + 1] - tables1d[1][ph_index]) / \
                      (tables2d[0][energy_index][ph_index + 1] - \
                       tables2d[0][energy_index][ph_index]) * \
                       (sigma - tables2d[0][energy_index][ph_index]);
    //cout<<"ek child is now "<<ek_ratio<<endl;
    child.energy *= ek_ratio;
    child.species = 3;
    child.charge = 0.0;
    child.mass = 0.0;
    //double costheta = calc_angular(parent.velocity);
    child.velocity = c_light_speed * normalize(parent.velocity);
    update_qed_factor(child);
    init_ptcls(child);
    init_ptcls(parent);
    /// 这个地方计算正负电子的速度应该根据动量计算，一会儿再修改：
    /// p = (1 - hv / cp ) P
    //这两项在boris　中更新
    double old_gamma = parent.lorentz_gamma;
    parent.energy = parent.energy - child.energy;
    parent.lorentz_gamma = parent.energy / mc2 + 1.0;
    
    parent.velocity = parent.velocity * old_gamma / parent.lorentz_gamma * (1.0 - child.energy / old_gamma / \
            parent.mass / c_light_speed / Mod(parent.velocity));
    if(child.energy < ekmin)
    {
        return 0;
    }
    return 1;

}
/***
double calc_angular(const __Vect3<double>& velocity)
{
    double ran1 = random_generator();
    double ran2 = random_generator();
    double ct = 2.0 * ran2 - 1.0;
    double beta = normalize(velocity) / c_light_speed;
    
***/

// 可能有2个表0,eta, 1chi, 2 f, 3 h(eta), 4 T(chi), 5, chiforT, 6, etaforh
// Pchi(eta, chi) Pf(chi, f), chi(eta, P)
// 返回值代表新的粒子个数

int bw_pair_create(const double& dt, __Ptcls& parent, __Ptcls& child0, __Ptcls& child1, \
        vector<vector<double> >& table, vector<vector<vector<double> > >& table_P)
{
    if(parent.rod_factor == 0)
    {
        //cout<<"initing rod"<<endl;
        init_ptcls(parent);
    }
    update_qed_factor(parent);
    if(parent.energy <= 2.0 * mc2)
    {
        //cout<<" energy too small "<<parent.energy<<" J, "<<parent.energy/MeV<<" MeV"<<endl;
        return 0;
    }
    double temp_qed = log10(parent.qed_factor);
    int chi_index_T = bs_search(temp_qed, table[5]);
    if(chi_index_T < 0 || chi_index_T >= table[5].size())
    {
        //cout<<"invalid chi_index_T "<<chi_index_T<<"\t"<<parent.qed_factor<<"\t field "\
            <<parent.bfield<<"\t"<<parent.efield<<"\t"<<Mod(parent.velocity) / c_light_speed<<endl;
        return 0;
    }
    /// 这里的系数先算成常数比较好
    parent.od_factor += 2.0 * PI * alpha_fine * e_mass * c_light_speed * c_light_speed /\
                        (lambda_compton * parent.energy) * parent.qed_factor *\
                        pow(10.0, table[4][chi_index_T]);
    if(parent.od_factor < parent.rod_factor)
    {
        return 0;
    }
    // Ｐ是chi 和eta的函数
    int chi_index = bs_search(temp_qed, table[1]);
    //cout<<" will generate pair with energy "<<parent.energy/MeV<<" MeV"<<endl;
    //cout<<" found chi_index "<<chi_index<<endl;
    if(chi_index < 0 ||chi_index >= table_P[1].size())
    {
        //cout<<" field strength too low "<<temp_qed<<endl;
        return 0;
    }
    int max_indP = table_P[1][chi_index].size() - 1;
    double random_P = random_generator() * table_P[1][chi_index][max_indP];
    int f_index = bs_search(random_P, table_P[1][chi_index]);
    //cout<<" findex is "<<f_index<<endl;
    while( f_index < 0 || f_index > table_P[1][chi_index].size())
    {
        random_P = random_generator() * table_P[1][chi_index][max_indP];
        f_index = bs_search(random_P, table_P[1][chi_index]);
    
    }
    double f_val = table[2][f_index];
    /// 通过从父粒子获取信息，因为很多东西共享
    //很明显ＩＤ这东西没法通用了
    child0 = parent;
    child1 = parent;
    child0.species = 0;
    child1.species = 2;
    child0.mass = e_mass;
    child1.mass = e_mass;
    child0.charge = - e_charge;
    child1.charge = e_charge;
    /// 从Pf中查找f给出正电子的能量。
    child0.energy = (parent.energy - 2.0 * mc2) * (1.0 - f_val);
    child1.energy = (parent.energy - 2.0 * mc2) * f_val;
    if(child0.energy < 0|| child1.energy < 0)
    {
        cout<<" pair creating error, negative energy!"<<endl;
    }
    child0.lorentz_gamma = child0.energy / mc2 + 1.0;
    child1.lorentz_gamma = child1.energy / mc2 + 1.0;
    child0.velocity = c_light_speed * sqrt(1.0 - pow(child0.lorentz_gamma, - 2.0)) \
                      * normalize(parent.velocity);
    child1.velocity = c_light_speed * sqrt(1.0 - pow(child1.lorentz_gamma, - 2.0)) \
                      * normalize(parent.velocity);
    init_ptcls(child0);
    init_ptcls(child1);
    update_qed_factor(child0);
    update_qed_factor(child1);
#ifdef __Debug
    if(std::isnan(child0.energy * child1.energy)) cout<<" nan electron or positron in bw pair "<<endl;
#endif
    //cout<<" positron vel "<<child1.velocity<<endl;
    return 2;
}


/// tables1d 0: ek,  1: x = omega / ek
/// tables2d: 0 sigma_br(ek_ele, x), 1, sigma_bh(ek_photon, y)
int bh_pair_create(const double& dt, const double& density, __Ptcls& parent, \
        __Ptcls& child0, __Ptcls& child1, vector<vector<double> >& tables1d, \
        vector<vector<vector<double> > >& tables2d)
{
    if(parent.energy < 2.0 * mc2)
    {
        return 0;
    }
    double ek_all_a = parent.energy;
    int ph_index = bs_search(parent.energy, tables1d[0]);
    //cout<<" ph f_index "<<ph_index<<endl;
    if(ph_index < 0)
    {
        return 0;
    }
    int max_ind = tables2d[0][0].size() - 1;
    //double sigma_interp = (log(tables2d[1][ph_index + 1][max_ind]) -  \
            log(tables2d[1][ph_index][max_ind])) / \
            (log(tables1d[1][ph_index + 1]) - log(tables1d[1][ph_index])) * \
            (log(parent.energy) - log(tables1d[1][ph_index])) + \
                          tables2d[1][ph_index][max_ind];
    /// 其实，最大截面与粒子的能量无关
    /***
    cout<<" ph f_index "<<ph_index<<endl;
    cout<<" parent ek "<<parent.energy<<" in table "<<tables1d[0][ph_index]<<endl;
    ***/
    double sigma_interp = tables2d[1][ph_index][max_ind];        
    double freq = Mod(parent.velocity) * density * sigma_interp * dt;
    double probability = 1.0 - exp(-freq);
    //cout<<" probability is "<<probability<<endl;
    if(random_generator() > probability)
    {
        return 0;
    }
    /// 决定产生粒子的能量,这个地方给出的是正电子能量占光子能量比
    double ek_all = parent.energy;
    double sigma_new = sigma_interp * random_generator();
    int ek_index = bs_search(sigma_new, tables2d[1][ph_index]);
    if(ek_index < 0)
    {
        return 0;
    }
    child0 = parent;
    child1 = parent;
    double ratio = (tables1d[1][ek_index + 1] - tables1d[1][ek_index]) / \
                   (tables2d[1][ph_index][ek_index + 1] - tables2d[1][ph_index][ek_index]) * \
                   (sigma_new - tables2d[1][ph_index][ek_index]) + tables1d[1][ek_index];
    ek_all = parent.energy - 2.0 * mc2;
    child1.energy = ek_all * ratio;
    child0.energy = ek_all - child1.energy;
    child0.mass = e_mass;
    child1.mass = e_mass;
    child0.charge = - e_charge;
    child1.charge = e_charge;
    child0.species = 0;
    child1.species = 2;
    init_ptcls(child0);
    init_ptcls(child1);
    child0.lorentz_gamma = child0.energy / mc2 + 1.0;
    child1.lorentz_gamma = child1.energy / mc2 + 1.0;
    child1.velocity = normalize(parent.velocity) * init_velocity(child1.lorentz_gamma);
    child0.velocity = normalize(parent.velocity) * init_velocity(child0.lorentz_gamma);
    return 2;
}



#endif
